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We consider the dynamics of finite-size disordered systems as defined by a master equation sat- 
isfying detailed balance. The master equation can be mapped onto a Schrodinger equation in 
configuration space, where the quantum Hamiltonian H has the generic form of an Anderson lo- 
calization tight-binding model. The largest relaxation time t eq governing the convergence towards 
Boltzmann equilibrium is determined by the lowest non- vanishing eigenvalue Ex = l/t eq of H (the 
, lowest eigenvalue being Eo — 0). So the relaxation time t eq can be computed without simulating 

■ the dynamics by any eigenvalue method able to compute the first excited energy Ex- Here we 

use the 'conjugate gradient' method to determine Ex in each disordered sample and present nu- 
i merical results on the statistics of the relaxation time t eq over the disordered samples of a given 

size for two models : (i) for the random walk in a self-affine potential of Hurst exponent H on a 
two-dimensional square of size L x L, we find the activated scaling lnt eq (L) ~ with ip = H as 

Q expected; (ii) for the dynamics of the Sherrington-Kirkpatrick spin-glass model of N spins, we find 
the growth \rit eq (N) ~ AT* with if) = 1/3 in agreement with most previous Monte-Carlo measures. 
t-H ' In addition, we find that the rescaled distribution of (lnt e<3 ) decays as e~" for large u with a tail 

exponent of order n ~ 1.36. We give a rare-event interpretation of this value, that points towards a 
' 1 ' . sample-to-sample fluctuation exponent of order %p W idth — 0.26 for the barrier. 

c ■ 

I. INTRODUCTION 

: 

' The non-equilibrium dynamics of disordered systems has been much studied both experimentally and theoretically 
(see for instance the reviews [l|, Q and references therein) . In numerical simulations, the main limitation is that the 
equilibrium time t eq (L) needed to converge towards equilibrium for a finite system of linear size L grows very rapidly 
^ \ with L. Within the droplet scaling theory proposed both for spin-glasses [H, HJ and for the directed polymer in a 
random medium [5j, the non-equilibrium dynamics is activated with barriers scaling as a power law B(L) ~ with 
some barrier exponent ip > that is independent of temperature and disorder strength. The equilibrium time t eq (L) 
O , then grows as 

CS,: lnteq(L) = B(L) ~ (1) 

This logarithmic scaling has been used to fit numerical data for disordered ferromagnets [6|,l7|,[8| and spin-glasses 

Sm- 
other authors, both for disordered ferromagnets and spin-glasses [13, Ell prefer a scenario corresponding to 

00 ■ logarithmic barriers B(L) ~ z(T, e)lnL, so that the equilibrium time t eq (L) scales as a power-law 

: t eq {L) = e B ^~L^ (2) 

Q\ \ where the exponent z(T, e) is non-universal and depends on the temperature T as well as on the disorder strength 
e. In the field of directed polymers or elastic lines in random media, the fit based the algebraic form of Eq. [2] used 
initially by many authors |l5[ has been now excluded by more recent work [l6|, [l?], EH , and has been interpreted as 
an artefact of an initial transient regime [13, EH • The reason why the debate between the two possibilities of Eqs [1] 
and [2] has remained controversial over the years for many interesting disordered models is that the equilibrium time 
t eq (L) grows numerically so rapidly with L that t eq (L) can be reached at the end of dynamical simulations only for 
rather small system sizes L < L max . For instance, in Monte-Carlo simulations of 2D or 3D random ferromagnets 

[1, 0, 1 EH, E2, El 

or spin-glasses Eol . EH E3, the maximal equilibrated size is usually only of order L max ~ 10 
lattice spacings. Even faster-than-the-clock Monte Carlo algorithms [20| . where each iteration leads to a movement, 
become inefficient because they face the 'futility' problem [2l|] : the number of different configurations visited during 
the simulation remains very small with respect to the accepted moves, i.e. the system visits over and over again 
the same configurations within a given valley before it is able to escape towards another valley. A recent proposal 
to improve significantly Monte Carlo simulations of disordered systems consists in introducing some renormalization 
ideas [22l ]. 

Taking into account these difficulties, a natural question is whether it could be possible to obtain informations on the 
equilibrium time t eq (L) without simulating the dynamics. In previous works [23l . |24| , we have proposed for instance 
to study the flow of some strong disorder renormalization procedure acting on the transitions rates of the master 
equation. However this approach is expected to become asymptotically exact only if the probability distribution of 
renormalized transitions rates flows towards an 'infinite disorder' fixed point, i.e. only for the activated scaling of Eq. 
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[TJ In the present paper, we test another strategy to compute t eq which is a priori valid for any dynamics defined by 
a master equation satisfying detailed balance : it is based on the computation of the first excited energy Ei of the 
quantum Hamiltonian H that can be associated to the master equation. This approach makes no assumption on the 
nature of the dynamics and is thus valid both for activated or non-activated dynamics (Eqs [T] or The mapping 
between continuous-time stochastic dynamics with detailed balance and quantum Schrodinger equations is of course 
very well-known and can be found in most textbooks on stochastic processes (see for instance (25|, [2(| H3] ) . However, 
since it is very often explained on special cases, either only in one-dimension, or only for continuous space, or only 
for Fokker-Planck equations, we stress here that this mapping is valid for any master equation satisfying detailed 
balance (see more details in section [TTJ) . In the field of disordered systems, this mapping has been very much used 
for one-dimensional models (see the review [28[ and references therein, as well as more recent works [H, HO, HH ) , but 
to the best of our knowledge, it has not been used in higher dimension, nor for many-body problems. In the field 
of many-body dynamics without disorder, this mapping has been already used as a numerical tool to measure very 
precisely the dynamical exponent z of the two dimensional Ising model at criticality 32]. 

The paper is organized as follows. In section|nl we recall how the master equation can be mapped onto a Schrodinger 
equation in configuration space, and describe how the equilibrium time t eq can be obtained from the associated 
quantum Hamiltonian. We then apply this method to two types of disordered models : section IIIII concerns the 
problem of a random walk in a two-dimensional self-affine potential, and section IIVI is devoted to the the dynamics 
of the Sherrington-Kirkpatrick spin-glass model. Our conclusions are summarized in section IVl 



II. QUANTUM HAMILTONIAN ASSOCIATED TO THE MASTER EQUATION 



A. Master Equation satisfying detailed balance 



In statistical physics, it is convenient to consider continuous-time stochastic dynamics defined by a master equation 
of the form 

= J2 p t (O w i c ' ^C)-P t (C) W out (C) (3) 
a 

that describes the the evolution of the probability Pt{C) to be in configuration C at time t. The notation W (C — » C) 
represents the transition rate per unit time from configuration C to C, and 

W out (C) = J2w(C^C) (4) 

c 

represents the total exit rate out of configuration C. Let us call U(C) the energy of configuration C. To ensure the 
convergence towards Boltzmann equilibrium at temperature T in any finite system 

U{C) 

P eq (C) = e -^- (5) 

where Z is the partition function 

Z = ^e-^ (6) 
c 

it is sufficient to impose the detailed-balance property 

U(C) ulc') 

e-TW(C^C')=e-TW{C'^C) (7) 



B. Mapping onto a Schrodinger equation in configuration space 



As is well known (see for instance [251 l2q . |27| ) the master equation operator can be transformed into a symmetric 
operator via the change of variable 

P f (C) = e-^ t (C) (8) 
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The function ipt{C) then satisfies an imaginary-time Schrodinger equation 

d ^l = -H M C) (9) 

where the quantum Hamiltonian has the generic form of an Anderson localization model in configuration space 

H = ^2e(C)\C ><C\+^2V{C,C')\C ><C'\ (10) 
c c,c 

The on-site energies read 

e(C) = W out (C) (11) 

whereas the hopping terms read 

V(C^C') = -e- Lm ^ m W{C'^C) (12) 

C. Specific choices for the detailed balance dynamics 

To have the detailed balance of Eq. it is convenient to rewrite the rates in the following form 

W (C _ C >) = 5 <c ,c<> e-^^e-s^) (13) 

where <5<c,c> means that the two configurations are related by an elementary dynamical move, and where S(C,C) is 
an arbitrary symmetric function : S(C,C) = S(C ,C). 

1. Simplest choice S(C,C') = 
To have the detailed balance property of Eq. [7j the simplest choice in Eq. [TBI corresponds to S(C,C) = 

W (C -> C) = 5 <c ,c> e 2T (14) 
Then the hopping terms of the quantum Hamiltonian are simply 

V(C,C) = -8 <c ,c>> (15) 

i.e. the non-vanishing hopping terms are not random, but take the same constant value (—1) as in usual Anderson 
localization tight binding models. The on-site energies are random and read 

e(C) = 2^ s <c,c> e ^ (16) 
c 

2. Metropolis choice 

In numerical simulations, one of the most frequent choice corresponds to the Metropolis transition rates 

r (u(c')-u(c)) ' 

W (C -> C) = 5 <e ,c>>mw l,e- ' r (17) 



In Eq. [13l this corresponds to the choice 



s^jun^m (18) 



In the quantum Hamiltonian, the hopping terms then read 
and the on-site energies are given by 



j. i- , \u{c')-u(c)\ 

V metr °P° hs (C,C') = -6 <c ,c>e ^ (19) 



metropolis 



(C) = J2 6 <W> min e (20) 
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D. Properties of the spectrum of the quantum Hamiltonian H 

Let us note E n the eigenvalues of H and \ip n > the associated normalized eigenvectors 

H\yj n > = E n \^ n > (21) 

J2\MC)\ 2 =1 (22) 
c 

The decomposition onto these eigenstates of the evolution operator e~ tH 

< C\e- tH \C Q >= £ e- B "Vn(C)C(Co) (23) 

a 

yields the following expansion for the conditional probability Pt (C\Cq) to be in configuration C at t if one starts from 
the configuration Co at time t = 

U(C) — UIC ) ,tj U(C)-U(C ) x — , „ . 

P t (C\C ) = e ^ < C\e~ tH \C >= e ^ ^ e'^^CW^) (24) 

n 

The quantum Hamiltonian H has special properties that come from its relation to the dynamical master equation : 
(i) the ground state energy is Eq = 0, and the corresponding eigenvector is given by 

p 2T 

MC) = (25) 



where Z is the partition function of Eq. [6l 

This corresponds to the convergence towards the Boltzmann equilibrium in Eq. [8] for any initial condition Co 

U{C) 

P t (C\C ) ~ e -" (C> - (Cn Vo(C)^(Co) = e -^— = Pe q (C) (26) 

t—* + oo Z 

(ii) the other energies E n > determine the relaxation towards equilibrium. In particular, the lowest non-vanishing 
energy Ei determines the largest relaxation time (1/Ei) of the system 

P t (C\C )-P eq (C) ~ e-^e- "^"^ MCWi(Co) (27) 

Since this largest relaxation time represents the 'equilibrium time', i.e. the characteristic time needed to converge 
towards equilibrium, we will use the following notation from now on 

te q = (28) 

The conclusion of this section is thus that the relaxation time t eq can be computed without simulating the dynamics 
by any eigenvalue method able to compute the first excited energy E\ of the quantum Hamiltonian H (where the 
ground state is given by Eq. [25] and has for eigenvalue Eq = 0). In the following subsection, we describe one of such 
methods called the 'conjugate gradient' method. 

E. Conjugate gradient method in each sample to compute E\ 

The 'conjugate gradient method' has been introduced as an iterative algorithm to find the minimum of functions 
of several variables with much better convergence properties than the 'steepest descent' method [H, [34| . It can be 
applied to find the ground state eigenvalue and the associated eigenvector by minimizing the corresponding Rayleigh 
quotient (35l.l36| 

R^<^ V > ( 29 ) 
< v\v > 
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The relation with the Lanczos method to solve large sparse eigenproblems is discussed in the chapters 9 and 10 of 
the book [34| . In the following, we slightly adapt the method described in (3f| [36[ concerning the ground state E 
to compute instead the first excited energy E\ : the only change is that the Rayleigh quotient has to be minimized 
within the space orthogonal to the ground state. 

In the remaining of this paper, we apply this method to various disordered models to obtain the probability 
distribution of the equilibrium time t eq {L) over the samples of a given size L. More precisely, since the appropriate 
variable is actually the equilibrium barrier defined as 



T e g — In t G q 



ln^i 



we will present numerical results for the probability distribution QL^eq) for various sizes L. 



(30) 



III. 



RANDOM WALK IN A TWO-DIMENSIONAL SELF-AFFINE POTENTIAL 





FIG. 1: (Color on line) Statistics of the equilibrium time t eq over the disordered samples of sizes I? for the random walk in 
a two-dimensional self-afline random potential of Hurst exponent H — 0.5 : (a) Probability distribution QL(T eq = lnt eg ) for 
L — 10, 20, 30, 40, 50, 60, 70, 80 (Inset : the corresponding distributions Q(u) of the rescaled variable u = (F eq — F eq (L))/ A(L) 
are shown in log scale for L — 10, 20, 30, 40); (b) the log-log plots of the disorder-average T eq {L) = In t eq (L) and of the width 
A(L) corresponds to the barrier exponent ijj = H — 0.5 (Eq. I34[) 

In this section, we apply the method of the previous section to the continuous-time random walk of a particle in a 
two-dimensional self-afhne quenched random potential of Hurst exponent H = 0.5. Since we have studied recently in 
[24j the very same model via some strong disorder renormalization procedure, we refer the reader to [24j and references 
therein for a detailed presentation of the model and of the numerical method to generate the random potential. Here 
we simply recall what is necessary for the present approach. 

We consider finite two-dimensional lattices of sizes Lx L. The continuous-time random walk in the random potential 
U(r) is defined by the master equation 

^±p. = P t ( f ') W(r'^ r) - P t (r) W out (r) (31) 

r ' 

where the transition rates are given by the Metropolis choice at temperature T (the numerical data presented below 
correspond to T = 1) 

W{?^r') = 5 < f if , > min(l,e-( c/ ( ? ')- c/ «)/ T ) (32) 

where the factor S < p_p '> means that the two positions are neighbors on the two-dimensional lattice. The random 
potential U(r) is self-affine with Hurst exponent H — 0.5 



[U(r)~U{r>)}\ c± \r-r'\ 2H (33) 

\f—r ' 1 
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On Fig. [T](a), we show the corresponding probability distribution <5i(r eg ) for various sizes 10 < L < 80 with a 
statistics of 36. 10 5 > n s (L) > 4.10 4 disordered samples. 

As shown by the log- log plots of Fig. [1] (b), we find that the disorder-averaged value T eq (L) and the width A(L) of 
the distribution Qz,(r e(? ) of the equilibrium barrier of Eq. 1301 involve the barrier exponent ip 



r e9 (L) 



L — >oo 



of value 



A(L) oc 

L — *-oo 



iP = H = 0.5 



(34) 



(35) 



These results are in agreement with scaling arguments on barriers [28l.l37| and with the strong disorder renormalization 
approach of [24| . 

IV. DYNAMICS OF THE SHERRINGTON-KIRKPATRICK SPIN-GLASS MODEL 



In(lnt) 
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FIG. 2: (Color on line) Statistics of the equilibrium time t eq over the disordered samples for the Sherrington-Kirkpatrick 
spin-glass model of N spins (2^ configurations) : (a) the log-log plot of the disorder-average F eq (L) as a function of TV for 
6 < N < 20 corresponds to the barrier exponent ip — 1/3 (Eq. I41[) (b) The rescaled probability distribution Q(u) of Eq. 1441 
shown here for 8 < N < 16, in log scale to see the tail of Eq. [45] : the tail exponent is of order r\ ~ 1.36. 

As an example of application to a many-body disordered system, we consider in this section the of the Sherrington- 
Kirkpatrick spin-glass model where a configuration C = {Si} of N spins Si = ±1 has for energy [38| 

U = - J v S i S i ( 36 ) 

l<i<j<N 

where the couplings are random quenched variables of zero mean J = and of variance J 2 = l/(N — 1). The 
Metropolis dynamics corresponds to the master equation of Eq. [3] in configuration space with the transition rates 

W (C - C) = 5 <c ,c '> min (l, e~^ c '^ U ^' T ) (37) 

where the factor <5<c.c '> means that the two configurations are related by a single spin flip. The data presented 
below correspond to the temperature T = 0.5 = T c /2. 

In the conjugate gradient described in section Til El one can start from a random trial vector to begin the iterative 
method that will converge to the first excited eigenvector. However, in the case of spin models where U is unchanged 
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if one flips all the spins Si — > —Si, one knows that the largest relaxation time will correspond to a global flip of all 
the spins. In terms of the quantum Hamiltonian associated to the dynamics discussed in section HTl this means that 
the ground state ipo of Eq. [25] is symmetric under a global flip of all the spins, whereas the first excited state ipi is 
anti-symmetric under a global flip of all the spins. As a consequence, we have taken as initial trial eigenvector for the 
conjugate gradient method the vector \v > defined as follows : denoting C pre f = {Sf } and C pre f = {— Sf } the 
two opposite configurations where the ground state tpo of Eq. [25] is maximal, one introduces the overlap between an 
arbitrary configuration C and C pre f 

N 

Q(C,C pref ) = J2s t Sr f (38) 

i=l 

and the vector 



v(C) = sgn(Q(C,C pref ))MC) (39) 

This vector is anti-symmetric under a global flip of all the spins and thus orthogonal to the ground state ipo- Moreover, 
it has already a small Rayleigh quotient (Eq. |29|) because within each valley where the sign of the overlap is fixed, it 
coincides up to a global sign with the ground state tpo of zero energy. So the non-zero value of the Rayleigh quotient 
of Eq. [25] only comes from configurations of nearly zero overlap Q. As a consequence it is a good starting point for 
the conjugate gradient method to converge rapidly towards the true first excited state ipi. 

We have studied systems of 6 < N < 20 spins (the space of configurations is of size 2^), with a statistics of 
10 7 > n s (N) > 1150 of independent disordered samples to compute the probability distribution QN^eq) of the 
largest barrier defined as 

T eq = lnt eq (40) 
As shown on Fig. EJa), we find that the disorder averaged equilibrium barrier scales as 



T eq {N) = lnt eq (N) oc with V^0.33 (41) 

This result is in agreement with theoretical predictions [H, H(| and with most previous numerical measures [4l], |42|, 
IHL IHL EE] ■ It is also interesting to consider the sample-to-sample fluctuation exponent ip W idth that governs the width 
of the probability distribution of the barrier 

A(N) ^ (TlqiN) ~ (T7 q (N)) 2 ) 1/2 cx JV*»««- (42) 



Although the disorder-average value has been much studied numerically [41|, [42J, |43|, |4J, |45j , the only measure of ip w idth 
we are aware of, is given by Bittner and Janke [45j 

VWtfc - 0.25 (43) 

With our numerical data limited to small sizes 6 < N < 20, we see already the expected behavior of the disorder- 
average of Eq. [41] as shown on Fig. [2] (a), but we are unfortunately not able to measure the exponent ip W idth of Eq. 
22] from the variance. 

However, as shown on Fig. HJb), the probability distribution <5jv(r e g) convergences rapidly towards a fixed rescaled 
distribution Q 

1 nl..- r.,-l%7(jv) 



g^'-AWn" 5 AW ) (44) 
We find that the rescaled distribution Q(u) presents at large argument the exponential decay 

lnQ(u) oc -v? (45) 

u — > + oo 

with a tail exponent of order 



■q ~ 1.36 



(46) 
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(on Fig. [Hb), a straight line would correspond to 77 = 1. Here we see a clear curvature indicating 77 > 1. The value 
of Eq. US] has been estimated via a three-parameters fit lnQ(u) ~ a — bu 71 for the data in the range u > 1). We are 
not aware of any theoretical prediction or any previous numerical measure of this tail exponent 77 to compare with. 
However, it should have an interpretation in terms of rare events. If the tail is due to rare samples that occur with 
some exponentially small probability of order e _ ( cst ) Ar ) but which present an anomalously large barrier of order N@ , 
the consistency equation for the powers of N in the exponentials read, using Eqs l44l and l45l 

(/3 - ipwidth) r\ = a (47) 

We may now consider the contribution of various types of rare events : 

(i) the anomalously ferromagnetic samples correspond to a = 2 ( with probability of order e -( cst ) N ; the N 2 random 
variables Jij will be all positive) and to j3 = 3/2 (instead of being finite, the local field hi = J2j Jij^j on s P m &i 
will be of order N 1 ' 2 ). If ipwidth = 1/3, the corresponding tail exponent is 77 = y- = 1.714 which we have measured 
elsewhere [46j for the case of the ferromagnetic Sherrington-Kirkpatrick model. Since here we measure a significantly 
different value, we believe that the rare events dominating the tail for the spin-glass Sherrington-Kirkpatrick model 
are not these ferromagnetic rare samples. 

(ii) in a typical sample, the distribution P(h) of the local field extends down to h = 0, with the linear behavior 
p(h) oc h as h — > 0) [U |H, |49|. However, with an exponentially small probability of order e~( cst ^ N , the N local 
fields of the sample will remain finite, i.e. bigger than some finite threshold hi > K , and the corresponding barrier 
will be anomalously large and of order N. These rare samples, that have an 'anormalously strong spin-glass order', 
in the sense that all local fields remain finite, thus correspond to the values a = 1 = /3 in Eq. 27] For instance, if 
^width = 1/3, the corresponding tail exponent reads rj = 5, whereas if ipwidth = 1/4, the corresponding tail exponent 
reads 77 = |. Our measure of Eq. [46] corresponds to 

ipwidth = 1 - - - 0.26 (48) 
77 

A tentative conclusion would thus be the following : at the small sizes that we can study, we cannot measure the 
width exponent ipwidth from the variance, but we can measure the tail exponent 77 that contains the information on 
ipwidth if one can properly identify the rare events that dominate the tail. In the spin-glass phase considered here, 
we believe that the rare events dominating the tail are the rare samples described in (ii) that have an 'anormalously 
strong spin-glass order', in the sense that all local fields remain finite, so that our measure of the tail exponent of Eq. 
[46] would point towards the value of Eq. [48] for the width exponent, which is actually very close to the value of Eq. 
Unmeasured by Bittner and Janke [HI from the variance for large sizes TV. These two indications suggest that ipwidth 
could actually be strictly smaller than the exponent ip = 1/3 governing the disorder-average value (Eq. 141]) . To the 
best of our knowledge, this question has never been raised for the barrier statistics, but it has been much discussed 
for the statistics of the ground state energy in the SK model (see [5(1 HH, HH HE H3| an d references therein) , where 
the sample-to-sample exponent W idth of the ground state energy (or the finite temperature free energy) is claimed to 
be either 9 W idth = 1/4 or Q width = 1/6, but is considered, in any case, to be smaller than the exponent 6 av = 1/3 that 
governs the correction to extensivity of the disorder average. A natural question is also whether the values ip = 1/3 
and ipwidth — 1/4 found in the statistics of the dynamical barrier are related to the exponents 9 av = 1/3 and B width 
that appear in the statistics of the ground state energy. 

V. CONCLUSION 

In this paper, we have proposed to use the mapping between any master equation satisfying detailed balance and a 
Schrodinger equation in configuration space to compute the largest relaxation time t eq of the dynamics via lowest non- 
vanishing eigenvalue E\ = l/t eq of the corresponding quantum Hamiltonian H (the lowest eigenvalue being Eq = 0). 
This method allows to study the largest relaxation time t eq without simulating the dynamics by any eigenvalue method 
able to compute the first excited energy E±. In the present paper, we have used the 'conjugate gradient' method 
(which is a simple iterative algorithm related to the Lanczos method) to study the statistics of the equilibrium time 
in two disordered systems : 

(i) for the random walk in a two-dimensional self-affine potential of Hurst exponent H 

(ii) for the dynamics of the Sherrington-Kirkpatrick spin-glass model of N spins. 

The size of vectors used in the 'conjugate gradient' method is the size Afc of the configuration space for the 
dynamics: for instance it is Mc = L 2 for the case (i) of a single particle on the two-dimensional square L x L and it is 
Afc — 2 N for (ii) containing N classical spins. We have shown here that the conjugate gradient method was sufficient 
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to measure the barrier exponents for these two models, but it is clear that it will not be sufficient for spin models in 
dimension d — 2 or d — 3 where the size of the configuration space grows as 2 L , and that it should be replaced by a 
quantum Monte-Carlo method to evaluate E± . For instance for the dynamics of the pure two dimensional Ising model 
at criticality studied in [32], the conjugate-gradient method used for squares L 2 of sizes L < 5 has been replaced for 
bigger sizes 5<L<15bya quantum Monte-Carlo method appropriate to compute excited states [55} . We thus 
hope that the same strategy will be useful in the future to compute the equilibrium time of disordered spin models in 
dimension d = 2. 
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